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Abstract 

The kinetic battery model is a popular model of the dynamic be¬ 
havior of a conventional battery, useful to predict or optimize 
the time until battery depletion. The model however lacks cer¬ 
tain obvious aspects of batteries in-the-wild, especially with re¬ 
spect to (i) the effects of random influences and (ii) the behavior 
when charging up to capacity bounds. 

This paper considers the kinetic battery model with bounded 
capacity in the context of piecewise constant yet random charg¬ 
ing and discharging. The resulting model enables the time- 
dependent evaluation of the risk of battery depletion. This is 
exemplified in a power dependability study of a nano satellite 
mission. 

1. Introduction 

The kinetic battery model (KiBaM) is a popular representation 
of the dynamic behavior of the state-of-charge (SoC) of a conven¬ 
tional rechargeable battery p4j , 25]. Given a constant load, it 
characterizes the battery SoC by two coupled differential equa¬ 
tions. Empirical evaluations show that this model provides a 
good approximation of the SoC across various battery types [20, 

p. 

The original KiBaM does not take capacity bounds into consid¬ 
erations, it can thus be interpreted as assuming infinite capacity. 
Reality is unfortunately different. When studying the KiBaM op¬ 
erating with capacity bounds, it becomes apparent that charging 
and discharging are not dual to each other. However, opposite 
to the discharging process, the charging process near capacity 
bounds has not received dedicated attention in the literature. 
That problem is attacked in the present paper. 

Furthermore, statistical results obtained by experimenting 
with real of-the-shelf batteries suggest considerable variances 
in actual performance [5|, likely rooted in manufacturing and 
wear differences. This observation asks for a stochastic re¬ 
interpretation of the classical KiBaM to take the statistically ob¬ 
served SoC spread into account on the model level, and this 
is what the present paper develops - in a setting with capacity 
bounds. It views the KiBaM as a transformer of the continu¬ 
ous probability distribution describing the SoC at any real time 
point, thereby also supporting uncertainty and noise in the load 
process. 

The approach presented not only enables the treatment of 
randomness with respect to the battery itself, but also makes 
it possible to determine the SoC distribution after a sequence 
of piecewise constant, yet random charge or discharge loads. 
We develop the approach in a setting with continuous random¬ 
ness so as to directly support normal (i.e. Gaussian), Weibull or 
exponential distributions. We apply the model to a case study 


inspired by a nano satellite currently orbiting the earth [17], 
for which we need to superpose it with a periodic deterministic 
charge load, representing the infeed from on-board solar panels. 

The resulting battery model can be viewed as a particular 
stochastic hybrid system [ 1 2 3| [6] 9]|28j, developed without 
discretizing time. It can (for instance) for any given real time 
point provide probabilistic guarantees about the battery never 
being depleted before. 

The genuine contributions of the paper are: (i) The inter¬ 
pretation of the KiBaM as a transformer of SoC distributions, 
(ii) developed without discretizing time, (iii) considering both 
charging and discharging in the context of capacity bounds, (iv) 
applied to a case study of a low earth orbiting satellite. 

Related work. Haverkort and Jongerden [ 21 ] review broad re¬ 
search on various battery models. They discuss stochastic bat¬ 
tery models [25|[7j which view the KiBaM for a given load as a 
stochastic process, unlike our (more accurate) view as a deter¬ 
ministic transformer of the randomized initial conditions of the 
battery. Furthermore, in this survey, the problem of charging 
bounds does not get dedicated attention. 

Battery capacity has been addressed only by Boker et al 
They considered a discretized, unbounded KiBaM together with 
a possibly non-deterministic and cyclic load process, synthesiz¬ 
ing initial capacity bounds to power the process safely. Hence, 
capacity is here understood as an over-dimensioned initial con¬ 
dition and not as a truly limiting charging bound. 

Random loads on a battery, generated by a continuous-time 
Markov chain, have been previously studied by Cloth et al 0 - 
Their setting cannot be easily extended by charging since they 
view the available and bound charge levels as two types of ac¬ 
cumulated reward in a reward-inhomogeneous continuous time 
Markov chain. 

An extension of the KiBaM to scheduling has been considered 
by Jongerden et al Jl9 ]. They compute optimal schedules for 
multiple batteries in a discretized setting with only discharging. 
This has been taken up and improved using techniques from the 
planning domain [12]. 

2. The Kinetic Battery Model 

The kinetic battery model is a mathematical characterization 
of the state of charge of a battery. It differs from an ideal energy 
source by incorporating the fact that not all the energy stored in 
a battery is available at all times. The stored energy is divided 
into two portions, the available charge and the bound charge. 
Only the available charge may be consumed immediately by a 
load supported by the battery and thereby behaves similar to an 
idealized source. As time passes, some of the bound charge is 
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Figure 1: Evolution of the state of charge as time passes (top) 
with the battery strained by a piecewise-constant load (bottom). 
The initially available charge decreases heavily due to the load 
400 but the restricted diffusion makes the bound charge de¬ 
crease only slowly up to time 10; after that the battery under¬ 
goes a mild recharge, followed by a strong recharge and a mild 
recharge at the end. At all times the bound charge approaches 
the available charge by a speed proportional to the difference of 
the two values. 



Figure 2: The two-well-model of the KiBaM with the available 
charge on the right (exposed directly to the load 7) connected 
to the bound charge on the left by a pipe of width p. 


converted into available charge and is thus free to be consumed. 
This effect is coined the recovery effect as the available charge 
recovers to some extend during periods of low discharge or no 
discharge at all. The recovery effect agrees with our experiences 
using batteries: For instance, when a cellphone switches off due 
to an apparently empty battery, it often can be switched back 
on after waiting a few minutes. The battery seems to have re¬ 
covered. This diffusion between available and bound charge can 
take place in either direction depending on the amount of both 
types of energy stored in the battery. Thus, while charging the 
battery, available charge is converted to bound charge. This be¬ 
havior is illustrated by Figure [l] 


Coupled differential equations. The KiBaM is often depicted 
as two wells holding liquid, the available charge and the bound 
charge well, interconnected by a pipe that represents the dif¬ 
fusion of the two types of charge, see Figure [2] Formally, the 
KiBaM is characterized by two coupled differential equations. 


m 

m 
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( 1 ) 

( 2 ) 


Here, the functions a(t) and b(t) describe the available and 
bound charge respectively, 7 is a load on the battery, p is the 
diffusion rate between both wells and c is the width of the avail¬ 


able charge well. Thus 1 — c is the width of the bound charge 
well. Intuitively a(t)/c and b(t)/(l — c) are the level of the fluid 
stored in the available charge well and the bound charge well, 
respectively. By defining 


c(l-c) 

we can rewrite 0 and 0 to 

a(t) = —7 + ck • b(t) — (1 - c)k • a(t) 

b(t) = (l-c)fc-a(t)-cfc-b(t). 

We will use this version of the KiBaM ODEs throughout this pa¬ 
per. 


Solving the equations. Using Laplace transforms the KiBaM 
ODE system can be solved, arriving at 

a t ,Aa 0 ,b 0 ) = q a (t)a 0 +r a (t)b 0 + s a (t)I 

b t ,i(a 0 ,b 0 ) = q b (t)a 0 + r b (t)b 0 + s b (t)I 

where a 0 and b 0 are the initial available and bound charge lev¬ 
els and the time-dependent coefficients of a 0 , b 0 and 7 in the 
equations can be expressed as 


q ait) = 

(1 -c)e~ kt + c 

q*(0 = 

—(1 — c)e~ kt + (1 — c) 

ra(t) = 

—c • e~ kt +c 

r b (t) = 

c • e~ kt +(l-c) 


(1 — c)(e~ kt — 1) 

Sait) = 

k 


(1 — c)(l — e~ kt ) 

s bit) = 

t • (1 c). 
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From the solution we can see that the KiBaM is affine in a 0 and 
b 0 (and also 7). Thus we can combine the two functions into 
one vector valued linear mapping 


K,, 
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’ q a (t) r a (t) s a (t) 
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When t is clear from context, we simplify the notation and 
drop the argument of these coefficients (without dropping the 
time-dependency). From now on we prefer semicolon notation 
[a; b] to denote column vectors. (All vectors appearing in this 
paper are column vectors.) Furthermore, whenever we com¬ 
pare two vectors, e.g., [a; b] < [a', b'], we interpret the order 
component-wise. 

Example 1 The function K can be used to approximate the final 
SoC in Figure^ffor k = 1/100, c = 1/2, and o denoting function 
composition) by 

[a; b] = K 4410 o ° Ki 5 ,- 6 oo ° ^ 30,-100 0 ^ 10,400 [5000; 5000] 

^ 1 ^ 44,100 0 ^ 15,-600 0 ^30,-100 [2002; 3998] 

^ 1^44,100 ° ^-15,-600 [4802; 4198] 

**Woo [10732; 7268], 
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and with the last step in more details (denoting e 100 by E), 
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On the last line, the first summands (with E) stand for the spread 
of the values when the recovery effect has not converged yet (as for 
t —> oo, E —> OJ. For c — 1/2 and 2 ero Zoad, the recovery effect 
makes the difference of the available charge and the bound charge 
converge to 0. However, for non-zero load I, it does not converge to 
0 but to I jk which explains the difference in the second summands. 

Battery Depletion. A standard application of KiBaM is to find 
out whether a task can be performed with a given initial state of 
charge without depleting the battery. A task is a pair (T, J) with 
T being the task execution time, and I representing the load, 
imposed for duration T. 

For an execution time T and a load I, we say that a battery 
with a SoC [a 0 ; b 0 ] > [0;0] powers a task (T,7) if 

K tJ [a 0 ;h 0 ] > [0;0] V 0 < t < T. 

Let us stress that the state of charge of the battery evolves 
in negative numbers in the same way as in positive numbers be¬ 
cause the differential equations do not have any explicit bounds. 
To rule out that the SoC of the battery goes into negative num¬ 
bers and returns back, we need a certain form of monotonicity. 

Example 1 (Cont.) We notice that neither the available nor the 
bound charge are monotonous in the standard sense. In Figure [I] 
the bound charge is not monotonous on [10,40], the available 
charge is not monotonous on [55,100]. However, for instance, on 
[40,55], available charge is the first to get above the value 9000 
(and never crosses the boundary again). 

Lemma 1 For any k,I g R, ^ e {<, >}, and 0 < t < T, 

Va 0 K, Vh 0 k : h t J (a 0 , h 0 ) di k —> a t J (a 0 , h 0 ) ^ k, 

Va 0 ii k, Vh 0 k : a t j(a 0 , h 0 ) :6 k —y a TI (a 0 , h 0 ) :ii k. 


10732 

7268 

-35 


Furthermore, the transformation law for random variables 
enables the construction of unknown density functions from 
known ones given the relation between the corresponding ran¬ 
dom variables. Formally, for every d -dimensional random vector 
X and every injective and continuously differentiable function 
g : R d —> R d , we can express the density function of Y := g(X) 
at value y in the range of g as 


/yOO = fx (g \y)) ■ det(j ? -i(y)) 


(4) 


where J f (x) denotes the Jacobian of a mapping / evaluated at 
x. Let us recall that the Jacobian of / is the matrix of the partial 
derivatives of the mapping /. 


Joint Density of the State of Charge. In order to consider the 
KiBaM as a stochastic object, it appears natural to consider the 
vector [ a 0 ;b 0 ;I ] as being random. This naturally reflects the 
situation where the initial state of the battery is subject to per¬ 
turbations due to manufacturing or wear variances, and so is 
its load. Therefore, we assume the initial SoC is expressed by 
random variables A 0 ,B 0 endowed with a joint probability den¬ 
sity function f 0 and that the load on the battery is expressed by 
a random variable I endowed with a probability density func¬ 
tion g. We assume that the random variables I and ( A 0 ,B 0 ) are 
independent. 

Example 2 A second running example addresses the random 
KiBaM. Instead of a single (Dirac) SoC, we now assume that 
the joint density f 0 of the charge is, say, uniform over the area 
[4,6.5] x [4,6.5] (below). 



Intuitively speaking, the first property states that the available 
charge is always the first to cross a bound, the second property 
states that when the available charge crosses a bound it never 
returns back (for a given load). 

As a direct consequence of Lemma [l] we can easily figure out 
whether the battery powers a task. 

Lemma 2 A battery with a SoC [a 0 ; b 0 ] > [0;0] powers a task 
C T, I) if and only if K t , [a 0 ;b 0 ] > [0;0]. 

3. The Basic Random KiBaM 

Some basic notions from probability theory are needed for 
the further development. Let f x and f XxY denote the density 
function of a random variable X and the joint density function 
of a pair of random variables (X, Y), respectively. 

The conditional density function f x \ Y of X given the occurrence 
of the value y of Y is defined as f x \ Y (.x\y) = Axy(^./)//y(/)- 
From this expression for f XxY , we obtain by marginalization the 
density function f x as 

r oo 

fxM= A|y(^l/)/y(y) dy. (3) 


We shall illustrate our findings how the SoC distribution evolves 
as the time passes on this particular example. 


Let ( A t ,B t ) denote the random variables expressing SoC after 
time T of constant (but random) load I. We are interested in the 
joint probability distribution of ( A T ,B T ) Thus, for a given time 
point T we want to establish the joint density function of the 
vector [. A t ;B t ] given by 


A Y 


B T 





A) 

B 0 

I 


(5) 


Expressing the joint density using direct application of the trans¬ 
formation law for random variables is not possible because the 
mapping is not invertible. However, using the fact that I and 
(A 0 ,B 0 ) are independent, it is still possible to use the transfor¬ 
mation law of random variables so that ultimately we arrive at 
an analytic characterization of the joint density of ( A T ,B T ). 

In the following we will derive the conditional density of 
( A t ,B t ) under the condition that I = i for some arbitrary but 
fixed value i. As g is known, we afterwards accommodate for 
the missing information about I via integration over the range 
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of I. Knowing that I = i eliminates one source of randomness, 
0 can be rewritten to 


Aj 

B t 
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" % r a 

_ Tb _ 


1 1 
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which is an invertible linear mapping and thus allows to express 
the joint density of ( A T ,B T ) in terms of the density of ( A 0 ,B 0 ) 
via the transformation law of random variables. 

Inverting this mapping K r>i results in 
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By a straightforward computation, the determinant of the Jaco¬ 
bian ofKl 1 is det/^-i = e kT . Note that it is constant in a, b , and 

T,i K r,i 

i, it only depends on T. Thus, using 0 we arrive at the joint 
density of ( A T ,B T ) conditioned by I = i 

/ r (a, b I i) —fj^ Ti [A 0 ;B 0 ] i. a i ft) 

=/ 0 (K-j[a; b]) • |e fct | 

where f KTi [A 0 -,B 0 ] denotes the joint density of the random vector 
K T}i [A 0 ;B 0 ]. According to (|3j) we arrive at the unconditional 
density over ( A T ,B T ) via integration on i: 


Lemma 4 A battery with SoC f 0 powers with probability p > 0 a 
task (T, g) if and only if 

r oo r oo 

f,fa, b) d b da > p. 

Jo Jo 

Example 2 (Cont.) Thanks to the lemma, it suffices to perform 
the integration on the densities displayed in the previous plots in 
this running example. The probability to power the tasks (20, g) 
is 1, for the task (60, g) it is just & 0.968. 

4. Bounded Recharging 

For the evaluation of the long-run state of a battery, a good 
understanding of the charging process is as essential as under¬ 
standing the discharging process. Both are well supported by 
the theory developed so far, and have occurred in our examples 
in the form of negative loads. What is not treated in the theory 
yet is a capacity bound of the battery which is a real constraint 
in most applications. To the best of our knowledge, charging in 
KiBaM while respecting its capacity restrictions has not been ad¬ 
dressed even in the deterministic case. This is what we are going 
to develop first, and then extend to the randomized setting. 

Let us assume that the battery has capacity d divided into 
capacity a max = c • d of the available charge well and capacity 
b m ax = (1 — c) • d of the bound charge well. 


Lemma 3 Let T be execution time and g be load density. For an 
initial SoC f 0 over ( A 0 ,B 0 ) and task (T, g), the joint distribution 
of (A t ,B t ) is absolutely continuous with density f T given by 


f T (x,y)= f 0 (K T J[x;y]ye kT ■ g(i) di. 


Example 2 (Cont.) We return to our example assuming the den¬ 
sity g of the load being uniform between [-0.1,0.1]. Based on 
the expression from Lemma [5] we can compute the SoC of the bat¬ 
tery after task (20, g), displayed on the left , and (60, g), displayed 
on the right. We arbitrarily chose the parameters c = 0.5 and 

p = 0.002. 




Risk of Depletion. Let us transfer the problem of battery de¬ 
pletion into the stochastic setting. We say that a density f Q is 
positive if for any a, b such that either a < 0 or b < 0 we have 
/ 0 (a, b) = 0. For an execution time T > 0 and a load density g, 
we say that the battery with positive SoC f Q powers with proba¬ 
bility p > 0 a task (T, g) if 

Pr [V0 < t <T : (A t ,B t ) > (0,0)] > p. 

Due to the monotonicity of KiBaM from Lemma [I] this is 
equivalent to observing the probability of being empty only at 
time T. From Lemma[3]we obtain the following. 


Charging at full available charge. A battery with empty avail¬ 
able charge can no longer support its task. Notably, a battery 
with full available charge does not behave dual to the behavior 
at empty available charge, just because a battery with full avail¬ 
able charge continues to operate, and we thus need to consider 
its further charging. In this discussion we do not consider crys¬ 
tallization effects reported for some batteries when overcharg- 
ing ]5l|26[. These phenomena do affect battery wear, not con¬ 
sidered in our modeling efforts thus far. 

When the available charge reaches its capacity a max = c • d 
and is still charged further by (high enough) incoming current, 
its value stays constant and only the bound charge increases due 
to diffusion. Hence, we know that a(t) = 0 and — = d and we 
can thus modify 0 to 


b(t)= P 



( 2 ’) 


Note that this equation describes the behavior of the battery 
at time t only if the load satisfies 




Since I is constant and the diffusion is decreasing over time, the 
inequality holds for all t, provided 


-I > b(0) = p 



K on 
1-cJ’ 


or equivalently if b( 0) > b tresh (7) where 

1 - c 

KeshU ) = frmax + 1 -• ( 6 ) 

P 

For a smaller initial bound charge, the standard differential 
equations 0 and 0 apply, until the available charge hits the 
capacity bound again. Finally, solving the ODE © yields 


Lemma 5 Let T be an execution time, I be a load, and b 0 be 
a bound charge such that b 0 > b tresh (I). A battery with a 
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SoC (a max , b 0 ) reaches after the task (T,7) the state of charge 
(a max , b T (b 0 )) where 

b T (b 0 ) = e- ckT b 0 +(l-e- ckT )-b max (7) 

and k again stands for p / (c • (1 — c)). 

We notice that the resulting bound charge b T (b 0 ) does not 
further depend on I, i.e. one cannot make the battery charge 
faster by increasing the charging current. Furthermore, for a 
fixed b 0 , the curve of t —> bfb 0 ) is a negative exponential start¬ 
ing from the point b 0 with the full capacity b max of the bound 
charge being its limit. Thus, Lemma [5] also reveals that the 
bound charge in finite time never gets full and there is no need 
to describe this situation separately. Finally, we denote anal¬ 
ogously by K r [a 0 ; b 0 ] = [a 0 ; £> r (b 0 )] the linear mapping de¬ 
scribing the behavior at the upper bound. 

Example 1 (Cont.) If we put an upper bound of 9000 to the bat¬ 
tery scenario from Figure [I] the battery ends up with a slightly 
smaller charge at time 100. 



Hitting the capacity bound. For a given constant load I, we 
have seen two types of behavior of the battery: (i) before it hits 
the available charge capacity and (ii) after it hits the capacity. 
The remaining question is when it hits the capacity limit. For a 
given initial state (a 0 , b Q ) and a load I, this amounts to solving 

b 0 ) ^maxj 

which in turn yields an equation 


where t is the largest solution of and t = T — t. 

The KiBaM evolution has one fundamental property we rely 
on heavily later: for any fixed time and load, it is monotonous 
with respect to starting SoC. 

Lemma 6 Let (a, b ) and (a 7 , b') be two SoCs. For every t > 0 
and for every I e R it holds that 

(a,b)<(a , ,b / ) =^> K°[a;b]<K°[a , ;b / ]. 


Example 1 (Cont.) By introducing the bounds in the first run¬ 
ning example, the computation of the final SoC changes only in 
the interval [40,55]. Here, instead o/K 15 _ 6 , we apply K- t _ e for 
the first t w 7.8 time units, followed by K 15 _ f . 


The computation of the time point t is problematic, as men¬ 
tioned before. An alternative to numerical approximation of the 
exact time point of crossing a bound is provided by the follow¬ 
ing observation: If time point t is fixed, we can check whether 
the available charge will exceed a max after t time units. In this 
case it is possible to determine the charging current I necessary 
for the available charge to hit a max exactly after t time units, 
i.e. solving a tJ (a 0 , b 0 ) = a max for I instead of t. Let us denote 
the solution of this equation by 7(a 0 , b 0 ) for a SoC (a 0 , b 0 ) and 
conclude that it can be computed by 


Tf 1 N Ta r a T . fl max 

I(a 0 ,b 0 ) = - a 0 - b 0 +-. 

S n S„ s„ 


Such a slower charging rate 7(a 0 , b 0 ) allows us to define a con¬ 
servative under-approximation of the exact KiBaM K° ; [a 0 ; b 0 ]. 
We replace the case K t o K fjJ [a 0 ; b 0 ] where the upper bound 
is reached by K r j( a0jbo )[a 0 ; b 0 ]. We will henceforth refer to this 
under-approximation by K^°. 


Lemma 7 For any SoC (a, b) we have [a; b ] < [a; b]. 

Together with Lemma|6] the above lemma tells us that, from the 
moment we used K^°first, we will never overshoot the exact 
behavior of the KiBaM with bounds K D . This fact is illustrated 
by the next example. 


u • e + v • t + w = 0 

where u = a 0 (1 — c) — b 0 c + (c + 1) • I/k, v = —Ic, and w = 
cd — a 0 c — b 0 c — (1 — c) • I/k. This can be solved as 



where W is the product log function. It has no closed form but 
can be arbitrarily numerically approximated [ 8 ]. 

Deterministic KiBaM with lower and upper bounds. All the 
previous building blocks allow us to express easily the SoC of 
a deterministic KiBaM after powering a given task (T,7) when 
considering battery bounds. We define it as the following func¬ 
tion: 


K T,I [ a o; b 0 ] 


Kr,j [a 0 ;b 0 ] 

° K i,i [ a oi b 0 ] 
[ 0 ; 0 ] 


v 


if a 0 > 0 and 
0 < a Tj/ (a 0 , b 0 ) < a max , 
if a 0 > 0 and 

a T,/( a 0? K)> a max? 

otherwise, i.e. if a 0 = 0 
or a rj (a 0 , b 0 ) < 0 


Example 1 (Cont.) We keep the upper bound of 9000 but instead 
of using the exact KiBaM behavior with bounds K n , we use the K^ n 
approximation. The load in the interval [40,55] is approximately 
-432.49 instead of -600 as the available charge would cross its 
bound. Instead it reaches 9000 exactly at t = 55. Note that from 
here on (i.e. in the interval (55,100]J the SoC is not exceeding 
the corresponding SoC from the previous figure. 



Random KiBaM with lower and upper bounds. We now turn 
our attention to the challenge of assuming that the random vari¬ 
ables ( A t ,B t ) evolve according to K° ; . We first observe that the 
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joint distribution of ( A T ,B T ) may not be absolutely continuous, 
because positive probability may accumulate in the point (0,0) 
where the battery is empty and on the line {(a max , b) | 0 < b < 
b max } where the available charge is full. Hence, we represent 
each ( A t ,B t ) by a triple {f tJ f t ,z t ) where 

• f t is the joint density describing the distribution in the 
“inner” area (0, a max ) x (0, b max ), 

• f t is the density over bound charge describing the distri¬ 
bution on the upper line {a max } x (0, b max ), and 

• z t G [0,1] is the probability of being empty 
such that for any measurable ACRxRwe have 

r a max r frmax 

Pr [CA t ,B t ) e A] = f t (a, b) l (a>b)eA db da 

Jo Jo 

+ I ft(b)l-(a max ,b)eA db +z t l(0,0)eA 

Jo 

where ± v denotes the indicator function of a condition <p. 

For an initial SoC ( f 0 ,f 0 ,z 0 ) and for a given task (T,g), our 
aim is to express the resulting SoC (/ r ,/ r ,z r ). 

First, we address the evolution inside the bounds by defin¬ 
ing f T . The new part is to describe how the battery moves 
from the upper bound to the area inside the bounds. We de¬ 
fine a vector valued function [b;i] K T>i [a max ; b~\ that maps 
(for a SoC with full available charge a max ) a pair of bound 
charge and load [b;i] to a SoC [a; b] in the usual KiBaM fash¬ 
ion. We denote the inverse of this mapping component-wise 
as [a;b] —> [^B(a, b); J(a, b)] where the functions ¥> (a, b) and 
J(a, b) can be explicitly expressed as 

J(a, b) = (a max e _fct - r b a - q b b)/ (v a s b - r b s a ) , 

■B(a, b) = -q 6 a + q a b + Cq 6 s a - q a s,,) • J(a, b ). 

The Jacobian determinant of this inverse map is easily derived 
to be 1/ (r a s fe — s a r b ) and is constant in the SoC and the load. 

The joint density has for any a < a max and b < b max the form 

/• oo 

fr( a >b) = I fo (k^J [a; b]) • |e fcT | • g(i) di (9) 

J —OO 

+ /oCB(a, b)) • |1/ (r a s b - s a r b ) | • g(J(a, b)). 

The first summand in 0 comes from the density f 0 of the 
inner area by the standard unbounded KiBaM. Ranging over all 
loads i, it integrates the density f 0 of such points (a^bj that 
satisfy K Tji [a f ;b f ] = [a; b], i.e. =K ~][a;b]. Lemmajl] 

again guarantees correctness that the bounds are not crossed in 
the meantime. The second summand comes from / 0 , due to dis¬ 
charging the battery down from the capacity bound as discussed 
above. Both summands are illustrated on the left hand side of 
Figure [3] 

After obtaining the density f T , we turn to f T , which is difficult 
to characterize, we are not aware of any closed-form expression 
achieving this. Hence we resort at this point to a conservative 
under-approximation of the battery charge. For a random SoC 
(A» B o) g iven by (foJo> z o) and a task (T,g) we define (A' T ,B' T ) 
given by {f^f^Zj). We define the densities for (0,0) < (a, b) < 



Figure 3: Illustration of the upper bound in random KiBaM. 
On the left, we depict the evolution towards a point (a, b) in 
the inner area of the state space. The figure shows different 
points from which the battery in the given time T reaches (for 
different loads) the point (a, b). Only one point corresponding 
to a unique load has available charge a max . On the right, we 
show three ways to reach the upper bound. 


C^max? ^max) 

f'M =Ma,b ) 

f* oo 

m= [ /o (b) • e cfcr • li>6 tresh( o 
J —OO 

+ fo(b) • lb<b tresh (i) • ^a T) i(a max ,b)>a max 

+ f fo (K- 1 [a max ; b] ) • e fct dt 1 • g(i) di 
Jo 

r 0 r 0 

z' T = I /r( a J b) db da 

J —oo J —oo 

where b denotes b^^b) = e ckT • b + (1 - e cfcr )(l - c)d. 

For expressing and z' T , we use the exact evolution given by 

t . Let us now closer discuss the density at the upper bound 
ich is an integral over all loads i. 

The first summand in the integral comes from the density f 0 of 
a point (a max , b) at the capacity bound such that K r [a max ; b] = 
[a ma x ‘,b]. This summand is taken into account only for such 
( a max> b) where the charging current covers the diffusion so that 
the battery evolves along the capacity bound as expressed by 
Lemma [ 5 ] Technically, we again apply the transformation law 
for random variables. 

Let us now address the case that the diffusion in a state 
( a max> b') at the upper bound is stronger than the charging cur¬ 
rent so that the available charge sinks in the beginning but 
before time T it again hits the upper capacity in some state 
( a max> b"). We are not able to express b' using b"; hence, we 
cannot “move” the density from b' to b". As apparent in the 
second summand, we thus underapproximate here the bound 
charge by assuming that the density stays in such state (a max , b'). 

The third summand in j' T comes from the density f Q of the 
inner area and is another under-approximation of bound charge. 
If available charge of the battery reaches the capacity bound 
before time T, we assume that until time T the SoC does not 
further evolve. In particular, the density goes to state (a max , b) 
from all states (a 0 , b 0 ) such that 

K t,i K; bo] = [“max; b] for some 0 <t<T. 

All three summands are illustrated in Figure [3] on the right. 
Let us finally state that ( A' T ,B' T ) is an under-approximation of 
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Figure 4: An MTP model on the left (with A depicted next to 
states) and its induced graph for T = 100 on the right. 


Lemma 8 For a random SoC ( A 0 ,B 0 ) given by (fo,f 0 ,z 0 ), a task 
(T, g) defining ( A T ,B T ) and ( A' T ,B' T ), and any fixed SoC (a, b), 

Pr([A r ;8 r ] > [a, b]) > Pr([A' T ;^] > [a, b]). 

Example 2 (Cont.) Based on Lemma [S] we can approximate the 
SoC of the random battery from our second running example for 
battery bounds [0,10]. We consider the same tasks , (20, g) on 
the left and (60, g) on the right 

).0622 


0.1 
0.08 
0.06 g 

0.04 

0.02 

0 

10308 


The bounded area of the joint density f T is depicted by the largest 
box. In the small box above we display the density f T at the ca¬ 
pacity limit. The numbers above and below are the probabilities of 
available charge being full and empty , respectively (the color below 
corresponds to the probability). 

5. Markov Task Process 

So far, we have only discussed execution of one task with fixed 
duration and random load. In this section, we give a discrete¬ 
time Markov model that randomly generates tasks that we call 
Markov task process (MTP). The formalism is closely inspired by 
stochastic task graph models [27] or data-flow formalisms such 
as SDF [22 ] or SADF [29]. In SDF, task durations are determin¬ 
istic, and thus directly supported in our framework. In SADF, 
durations are in general governed by discrete probability distri¬ 
butions, which can be translated into our framework at the price 
of a larger state space. 



(a 0 , b 0 ) of the battery and an initial task s 0 e S are chosen in¬ 
dependently at random according to {f 0 ,f 0 ,z 0 ), and n, respec¬ 
tively. Then, the load i 0 in task s 0 is picked randomly according 
to g(s 0 ). After the battery is strained by the load i 0 for A(s 0 ) 
time units, the process moves into a random successor task s 1 
(where any s 1 is chosen with probability P(s 0 ,s 1 )). Here, the 
load f is randomly chosen and so on. 

Formally, M and (/ 0 ,/ 0 ,2 0 ) induce a probability measure Pr 
over samples of the form co = [(a 0 , b 0 ); (s 0 , i 0 X 5 i> h )' * * ] where 
the first component is the initial SoC of the battery and the sec¬ 
ond component describes an infinite execution of M. Here, each 
Sj is the j-th task and ij is the load that is put on the battery for 
A (s_y) time units while performing Sj. For a given T e [R> 0 , the 
SoC of the battery at time T is expressed by random variables 
A t ,B t that are for any co = [(a 0 , b 0 ); (s 0 , i 0 )( s iHi)" * ] defined 
as 


Aj- (to) 

B r (a>) 


= K°. oK° . o • • • o K? . 

t ,l n f n-l4n-l 1 04c 


^0 


where each tj stands for A(s ; ), and n is the minimal number 
such that the n-th task is not finished before T, i.e. t n > t' 
where t' := T - tj. 


Definition 2 We say that a battery with SoC (/ 0 ,/ 0 ,z 0 ) powers 
with probability p > 0 a system M for time T if 

Pr [A t > 0] > p. 


In order to (under-) approximate the probability that M is 
powered for a given time, we need to symbolically express the 
distribution over ( A T ,B T ). We present an algorithm that builds 
upon the previous results. 


Expressing the distribution of ( A T ,B T ). Let us fix an input 
MTP M = (S, P, 7i, A, g), distribution over SoC (f 0 ,fo,z 0 ), and 
time T > 0. We consider the joint distribution of SoC and the 
MTP Intuitively, we split the distribution of SoC into subdistri¬ 
butions and move them along the paths of M according to the 
probabilistic branching of the MTP We notice that we do not 
need to explore all exponentially many paths; when two paths 
visit the same state at the same moment, we can again merge 
the two subdistributions. This process is formalized by the fol¬ 
lowing graph and a procedure how to propagate the distribution 
through the graph. 

For a given MTP M we define a directed acyclic graph (V,E) 
over V = S x ({0, 1,..., [T\] U {T}) such that there is an edge 
from a vertex (s, t) to a vertex (s', t') if P(s,s') > 0, t < t', and 
t' = min{t + A(s), T}. Further, let (V',E') be the graph obtained 
from (V, E) by removing vertices that are not reachable from any 
(s,0) with n (s) > 0 (see Figure|4|. 

1. We label each vertex of the form (s, 0) where 7i(s) > 0 by 
a subdistribution (f 0 • n(s),f 0 • tt(s),z 0 • tt(s)). 


Definition 1 A Markov task process (MTP) is a tuple M = 
(S,P, 7T, A, g) where S is a finite set of tasks, P : S x S —> [0,1] 
is a probability transition matrix , n is an initial probability dis¬ 
tribution over S, A : S —> N assigns to each task an integer time 
duration, and g assigns to each task a probability density function 
of the load. 


2. We repeat the following steps as long as possible. 

(a) For each vertex (s, t) labeled by ( f,f,z ), we obtain 
( f',f',z') by Lemmalijfor a task (t' - t,g(s)) where 
t' = min{t + A(s), T Then we label each outgoing 
edge from (s, t) to some (s', t') by 

(/' • P(s,s'), S' • P(s,s'), z' • Pis,s')} . 

(b) For each vertex (s, t) such that all its ingoing edges 
are labeled by (f 1 ,f 1 ,z 1 ),...,(f k ,f k ,z k ), we label 


An example of a MTP is depicted in Figure [4] Intuitively, a 
Markov task process M together with an initial distribution over 
SoC given by (fo,f 0 ,z 0 ) behaves as follows. First, an initial SoC 
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0 , t) by 

(f 1 + ---+f k ,P + ---+f k , z l + --- + z k ). 

Finally, let all vertices of the form (s, T) e V' be labeled by 
{f 1 ,f 1 ,z 1 ),...,{f n ,f n ,z n ). The output distribution (f T ,f T ,z T ) 
is equal to 

(z 1 + •••+/", f 1 + •••+/", Z 1 + • • • +Z n ) . 

We naturally arrive at the following theorem. 

Theorem 1 A battery with SoC (f 0 ,fo> z o) powers a system M for 
time T with probability at least 1 — z T . 

6. The Random KiBaM in Practice 

In this section, we apply the results established in the previ¬ 
ous sections in a concrete scenario. The problem is inspired by 
experiments currently being carried out with an earth orbiting 
nano satellite, the GOMX-1 ED- 

Satellite. GOMX-1 [17 ] is a Danish two-unit CubeSat mission 
launched in November 2013 to perform research and experi¬ 
mentation in space related to Software Defined Radio (SDR) 
with emphasis on receiving ADS-B signal from commercial air¬ 
craft over oceanic areas. As a secondary payload the satellite 
flies a NanoCam C1U color camera for earth observation exper¬ 
imentation. Five sides are covered with NanoPower PI 10 solar 
panels, the power system NanoPower P31u holds a 7.4V Li-Ion 
battery of capacity 5000 mAh. GOMX-1 uses a radio amateur 
frequency for transmitting telemetry data, making it possible to 
receive the satellite data with low-cost infrastructure anywhere 
on earth. The mission is developed in collaboration between 
GomSpace ApS, DSE Airport Solutions and Aalborg University, 
financially supported by the Danish National Advanced Technol¬ 
ogy Foundation. The empirical studies carried out with GOMX-1 
serve as a source for parameter values and motivate the sce¬ 
nario described in the sequel. We concretely use the following 
data collected from extensive in-flight telemetry logs. 

• One orbit takes 99 minutes and is nearly polar; 

• The battery capacity is d = 5000 mAh; 

• During 4 to 7 out of on average 15 orbits per day, com¬ 
munication with the base station takes place. The load 
induced by communication is roughly 400 mA. The length 
of the communication depends on the distance of the pass 
of the satellite to the base station and varies between 5 
and 15 minutes; 

• In each communication, the satellite can receive instruc¬ 
tions on what activities to perform next. This influences 
the subsequent background load. Three levels of back¬ 
ground load dominate the logs, with average loads at 250 
mA, 190 mA, and 90 mA. These background loads sub¬ 
sume the power needed for operating the respective ac¬ 
tivities, together with basic tasks such as sending beacons 
every 10 seconds; 

• Charging happens periodically, and spans around 2/3rd of 
the orbiting time. Average charge power is 400 mA; 

The above empirical observations determine the base line of our 
modeling efforts. Still the case study described below is a syn¬ 
thetic case. We make the following assumptions: 


• We assume constant battery temperature. The factual 
temperature of the orbiting battery oscillates between - 
8 and 25 degree Celsius on its outside. There is the 
(currently unused) on-board option to heat the battery to 
nearly constant temperature. Using an on-off controller, 
this would lead to another likely nearly periodic load on 
the battery, well in the scope of what our model supports. 

• A constant charge from the solar panels is assumed when 
exposed to the sun. The factual observed charge slowly 
decays. This is likely caused by the fact that solar panels 
operate better at lower temperature (opposite to batter¬ 
ies), but heat up quickly when coming out of eclipse. 

• We assume a strictly periodic charging behavior. The fac¬ 
tual charging follows a more complicated pattern deter¬ 
mined by the relative position of sun, earth and satellite. 
There is no fundamental obstacle to calculate and incor¬ 
porate that pattern. 

• We assume a uniform initial charge between 70% and 90% 
of full capacity with identical bound and available charge. 
Since the satellite needs to be switched off for transporta¬ 
tion into space, assuming an equilibrated battery is valid. 
Being a single experiment, the GOMX-1 had a particular 
initial charge (though unknown). The charge of the orbit¬ 
ing battery can only be observed indirectly, by the voltage 
sustained. 

• We assume that the relative distance to a base station is a 
random quantity, and thus interpret several of the above 
statistics probabilistically. In reality, the position of the 
base station for GOMX-1 is at a particular fixed location 
(Aalborg, Denmark). Our approach can either be viewed 
as a kind of probabilistic abstraction of the relative satel¬ 
lite position and uncertainty of signal transmission, or it 
can be seen as reflecting that base stations are scattered 
around the planet. This especially would be a realistic 
in scenarios where satellite-to-satellite communication is 
used. 

• We assume that the satellite has no protection against bat¬ 
tery depletion. In reality, the satellite has 2 levels of soft¬ 
ware protection, activated at voltage levels 7.2 and 6.5, re¬ 
spectively, backed up by a hardware protection activated 
at 6 V In these protection modes, various non-mission- 
critical functionality is switched-off. Despite omitting such 
power-saving modes, we still obtain conservative guaran¬ 
tees on the probability that the battery powers the satel¬ 
lite. 

Satellite model. According to the above discussion, the load 
on the satellite is the superposition of two piecewise constant 
loads. 

• A probabilistic load reflecting the different operation 
modes, modeled by a Markov task process Mas depicted 
in Figure [5] 

• A strictly periodic charge load alternating between 66 
minutes at —400 mA, and the remaining 33 minutes at 
0 mA. 

Actually it is after 364 days, as this is in the middle of the 
charging phase. After 365 days the satellite is in eclipse and no 
density is exibited in the upper diagram. 
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Figure 6: SoC of the satellite after 1 yearj^vith different sizes of the battery The leftmost SoC is with the original battery capacity 
5000 mAh. In each further plot, the battery capacity is halved, i.e. 2500 mAh, 1250 mAh, and 625 mAh. Note that all the densities are 
depicted on the logarithmic scale (numbers in the legend stand for the order of magnitude). We observe that only the smallest battery 
does not give sufficient guarantees. Its probability of depletion after 1 year is 0.0365; the probability gets down to 1.7 • 10 -10 already 
for the 1250 mAh battery. The smaller the battery, the more crucial is the distinction of available and bound charge as a larger area of 
the plots is filled with non-trivial density. 



5 min. 400 mA 


Figure 5: Markov task process of the load on the satellite. All 
load distributions are normal with mean depicted next to the 
states and with standard deviation 5. This load is superposed 
with a strictly periodic load modelling charge by solar power 
infeed. 


One can easily express the charging load as another independent 
Markov task process (where all probabilities are 1) and consider 
the sum load generated by these two processes in parallel (meth¬ 
ods in Section [5] adapt straightforwardly to this setting). 

The KiBaM in our model has following parameters: 


• the ratio of the available charge c — 1/2 (artificially cho¬ 
sen value as parameters fitted by experiments on similar 
batteries strongly vary [[31,20 ]); 


the diffusion rate p = 0.0006 per minute (we decreased 
the value reported by experiments [20] by a factor of 4 
because of the low average temperature in orbit, 3.5°C, 
and the influence of the Arrhenius equation J23]]). 


Computational Aspects. We implemented the continuous so¬ 
lution developed in the previous sections in a high-level compu¬ 
tational language Octave. This showed up to be practical only up 
to sequences of a handful of tasks. Therefore, we implemented a 
solution over a discretized abstraction of the stochastic process 
induced by the MTP and the battery. By fixing the number of 
discretization steps K e N which yields the discretization step 
6 = | we obtain battery states 

• {(n, m) \ 0 < n < K, 0 < m < K} in the inner space where 
each (n,m) represents the adjacent rectangle of higher 


charge [nd, (n + 1)6) x [m<5, (m + 1)6), 

• states {( K,m ) | 0 < m < K} on the capacity boundary 
where each (IC, m) represents the adjacent line of higher 
bound charge {|} x [md, (m + 1)6), and 

• state (0,0) for the rest, {(a, b) \ a < 6 or b < 6}. 

We always represent higher charge by lower charge (i.e. 
under-approximate SoC) since we are interested in guarantees 
on probabilities that the battery powers the MTP for a given time 
horizon. Similarly, we replace load distributions by discrete dis¬ 
tributions where each point represents an adjacent left interval 
(i.e., we over-approximate the load). The continuous methods 
of Lemma [8] are easily adapted to this discrete setting, basically 
replacing integrals by finite sums. 

This methods gives us an underapproximation of the proba¬ 
bility that the battery powers the satellite. We do not have any 
prior error bound, but one can make the results arbitrarily pre¬ 
cise by increasing K, at the price of quadratic cost increase. 

Our implementation is done in C++, we used K = 1200, 600, 
300 and 150 for the experiments with the batteries of capac¬ 
ity 5000 mAh, 2500 mAh, 1250 mAh and 625 mAh, respec¬ 
tively to guarantee equal relative precision. All the experiments 
have been performed on a machine equipped with an Intel Core 
i5-2520M CPU @ 2.50GHz and 4GB RAM. All values occuring 
are represented and calculated with standard IEEE 754 double¬ 
precision binary floating-point format except for the values re¬ 
lated to the battery being depleted where we use arbitrary preci¬ 
sion arithmetic (as to this number, we keep adding values from 
the inner area that are of much lower order of magnitude). 

Model evaluation. We performed various experiments with 
this model, to explore the random KiBaM technology. We here 
report on four distinct evaluations, demonstrating that valuable 
insight into the model can be obtained. 

1. The 5000 mAh battery in the real satellite is known to be 
over-dimensioned. Our aim was to find out how much. 
Hence, we performed a sequence of experiments, decreas¬ 
ing the size of the battery exponentially. The results are 
displayed and explained in Figure [6] We found out that 
1 /4 of the capacity still provides sufficient guarantees to 
power the satellite for 1 year while 1/8 of the capacity, 
625 mAh, does not. 
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Figure 7: Load noise. The 1 year run of the 1250 mAh battery 
with Dirac loads on the left and with noisy loads on the right. 
We used Gaussian noise with standard deviation 5. 



Figure 8 : Number of solar panels. The full 5000 mAh battery 
with 9 solar panels on the left and 6 solar panels on the right. 
This shows that for the current load, a 1 unit cube design with 
solar panels on 6 or less sides is not possible. 


2. We compared our results with a simple linear battery 
model of the same capacity. This linear model is not un¬ 
common in the satellite domain, it has for instance been 
used in the Envisat and CryoSat missions (15]. We obtain 
the following probabilities for battery depletion: 


capacity 

linear battery model 

KiBaM 

5000 mAh 

1 .86-10“ 84 

1.7-10“ 63 

625 mAh 

2.94-10“ 8 

0.0365 


The linear model turns out to be surprisingly (and likely 
unjustifiably) optimistic, especially for the 625 mAh bat¬ 
tery. 

3. We (computationally) simplified the two experiments 
above by assuming Dirac loads. To analyze the effect of 
the white noise, we compared the Dirac loads with the 
noisy loads, explained earlier, on the 625 mAh battery. As 
expected, the noise (a) smoothes out the distribution a lit¬ 
tle and (b) pushes a bit more of the distribution to full and 
empty states, see Figure [7] 

4. Our reference satellite is a two-unit satellite, i.e. is built 
from two cubes, each 10 cm per side. In the current de¬ 
sign, 9 of the 10 external sides are covered by solar pan¬ 
els, the remaining one is used for both radio antenna and 
camera. We thus analyzed whether a one-unit design with 
only 5 solar panels is possible. The answer is negative, the 
system runs out of energy rapidly with high probability. 
Figure [8] displays that even for 6 panels the charge level is 
highly insufficient to sustain the load. 


7. Alternative Approaches 

The results reported above are obtained from a discretized 
abstraction of the stochastic process induced by the MTP and the 
battery, solved numerically and with high-precision arithmetic 
where needed. 

One could instead consider estimating the probability z t of 
the battery depletion using ordinary simulation techniques [ 16]. 
Considering a battery of capacity 5000 mAh, this would mean 
that about 10 63 simulations traces are needed on average to ob¬ 
serve the rare event of a depleted battery at least once. This 
seems prohibitive, also if resorting to massively parallel simu¬ 
lation, which may reduce the exponent by a small constant at 
most. A possible way out of this might lie in the use of rare 
event simulation techniques to speed up simulation [30]. 

The behaviour of KiBaM with capacity bounds can be ex¬ 
pressed as a relatively simple hybrid automaton model [18]. 
Similarly, the random KiBaM with capacity bounds can be re¬ 
garded as an instance of a stochastic hybrid system (SHS) 
3][6|[9][28t. This observation opens some further evaluation av¬ 
enues, since there are multiple tools available publicly for check¬ 
ing reachability properties of SHS. In particular, Faust 2 0. 
Si Sat © and ProHVer (32 13j appear adequate at first sight. 
Our experiments with Faust 2 however were unsuccessful, ba¬ 
sically due to a model mismatch: The tool thus far assumes 
stochasticity in all dimensions, because it operates on stochastic 
kernels, while our model is non-stochastic in the bound charge 
dimension. With SiSat, we so far failed to encode the MTP (or 
its effect) into an input accepted by the tool. The MTP can be 
considered as a compact description of an otherwise intricate 
semi-Markov process running on a discrete time line. This is in 
principle supported by SiSat, yet we effectively failed to pro¬ 
vide a compact encoding. Our ProHVer experiments failed for 
a different reason, namely the sheer size of the problem. All the 
above tools have not been optimized for dealing with very low 
probabilities as they appear in the satellite case. 


8. Conclusion 

Inspired by the needs of an earth-orbiting satellite mission, 
we extended in this paper the theory of kinetic battery mod¬ 
els in two independent dimensions. First, we addressed battery 
charging up to full capacity. Second, we extended the theory of 
the KiBaM differential equations to a stochastic setting. We pro¬ 
vided a symbolic solution for random initial SoC and a sequence 
of piecewise-constant random loads. 

These sequences can be generated by a stochastic process rep¬ 
resenting an abstract and averaged behavioral model of a nano 
satellite operating in earth orbit, superposed with a determin¬ 
istic representation of the solar infeed in orbit. We illustrated 
the approach by several experiments performed on the model, 
especially varying the size of the battery, but also the number of 
solar panels. 

ESA is running a large educational program JlO] for launch¬ 
ing missions akin to GOMX-1. The satellites are designed by 
student teams, have the form of standardized 1 unit cube with 
maximum mass of 1 kg, and target mission times of up to four 
years. The random KiBaM presented here is of obvious high rele¬ 
vance for any participating team. It can help quantify the risk of 
premature depletion for the various battery dimensions at hand, 
and thereby enable an optimal use of the available weight and 
space budget. Our experiments show that using the simpler lin¬ 
ear battery model instead is far too optimistic in this respect. 

For a fixed setup, one can also use the technology offered 
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by us for optimal task scheduling: In the same way as we can 
follow a single SoC distribution, we can also branch into several 
distributions and determine which of them is best according to 
some metric. Taking inspiration from [ 3lJ, this can be combined 
with statistical model checking so as to find the optimal task 
schedule of a given set of tasks. 

Several extensions can and should be integrated in the model. 
Among them, temperature dependencies are of particular inter¬ 
est. A temperature change has namely opposing physical ef¬ 
fects in solar panels and in the battery, having intriguing con¬ 
sequences such as piecewise exponential decay in the charging 
process. An extension that is particularly important for long last¬ 
ing missions, is incorporating a model of battery wearout. So far 
we assume the battery capacity to be constant along the mission 
time. Notably, our contribution is the first to consider capacity 
bounds in operation at all, as far as we are aware. 
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